////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 12 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "Burgers" solves the 1D Burgrs equation using the Lax and Lax-Wendroff methods.           //
//  Periodic boundary condition is used. The system length L is set to one.                               //
//  Initial wave shape is Gaussian pulse.                                                                 //                         //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

double tau=0.001,h,L=1,sigma=0.05,x;

int N=400,i,n,T=10000;


ofstream file0 ("initial profile n=0.plt");  //output file for the profile at timestep n=0.
ofstream file1 ("profile Lax n=4500.plt");    //output file for the profile at timestep n=100.
ofstream file2 ("profile Lax n=5000.plt");    //output file for the profile at timestep n=500.
ofstream file3 ("profile Lax n=6000.plt");    //output file for the profile at timestep n=500.
ofstream file4 ("profile LW n=4500.plt");     //output file for the profile at timestep n=500.
ofstream file5 ("profile LW n=5000.plt");    //output file for the profile at timestep n=500.
ofstream file6 ("profile LW n=6000.plt");    //output file for the profile at timestep n=500.



vector<double> u(N+1,0),unew(N+1,0),uL(N+1,0),uLnew(N+1,0),uLW(N+1,0),uLWnew(N+1,0),J(N+1,0);

 // Arrays "u" and "unew" store the current and updated values
 // of the scaler function u(x,t) at grid points. "L" stands for Lax and "LW" stands for Lax-Wendroff.

 h=L/N; cout<<"h= "<<h<<endl; cout<<"tau= "<<tau<<endl; //getch();


for(i=0;i<=N;i++){

x=-0.5*L+i*h; u[i]=exp(-0.5*(x/sigma)*(x/sigma));

uL[i]=exp(-0.5*(x/sigma)*(x/sigma)); uLW[i]=exp(-0.5*(x/sigma)*(x/sigma));

file0<<(-0.5*L+i*h)<<"  "<<u[i]<<endl; }

//storing the initial condition data in the data file.


for (n=1;n<=T;n++){ //cout<<"n= "<<n<<endl; getch();



/*
// -------------------       FTCS method -------------------------------------

for(i=1;i<N;i++){ unew[i]=u[i]-(tau/(4*h))*( u[i+1]*u[i+1]-u[i-1]*u[i-1] ); }

unew[0]=u[0]-(tau/(4*h))*( u[1]*u[1]-u[N]*u[N] ); // periodic condition has been used: u[-1] is replaced with u[N]

unew[N]=u[N]-(tau/(4*h))*( u[0]*u[0]-u[N-1]*u[N-1] ); // periodic condition has been used: u[N+1] is replaced with u[0]


//if (n==10)  for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.

//if (n==20) for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.
*/





//-----------------------------     Lax Method      ------------------------------------

for(i=1;i<N;i++){ uLnew[i]=0.5*(uL[i+1]+uL[i-1])-(tau/(4*h))*( uL[i+1]*uL[i+1]-uL[i-1]*uL[i-1]); }

uLnew[0]=0.5*(uL[1]+uL[N])-(tau/(2*h))*( uL[1]*uL[1]-uL[N]*uL[N]);// periodic condition has been used: uL[-1] is replaced with uL[N]

uLnew[N]=0.5*(uL[0]+uL[N-1])-(tau/(2*h))*( uL[0]*uL[0]-uL[N-1]*uL[N-1]);// periodic condition has been used: uL[N+1] is replaced with uL[0]


if (n==4500)   for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.

if (n==5000)   for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.

if (n==6000)  for(i=0;i<=N;i++){ file3<<(-0.5*L+i*h)<<"  "<<uLnew[i]<<endl; } //storing in the data file.



//-----------------------------     Lax-Wendroff Method      ------------------------------------

for(i=0;i<=N;i++){ J[i]=0.5*uLW[i]*uLW[i]; }


uLWnew[0]=uLW[0]-(tau/(2*h))*(J[1]-J[N]) + (tau*tau/(4*h*h))*( (uLW[1]+uLW[0])*(J[1]-J[0])-(uLW[0]+uLW[N])*(J[0]-J[N]) );


for(i=1;i<N;i++){

uLWnew[i]=uLW[i]-(tau/(2*h))*(J[i+1]-J[i-1]) + (tau*tau/(4*h*h))*( (uLW[i+1]+uLW[i])*(J[i+1]-J[i])-(uLW[i]+uLW[i-1])*(J[i]-J[i-1]) );

}


uLWnew[N]=uLW[N]-(tau/(2*h))*(J[0]-J[N-1]) +  (tau*tau/(4*h*h))*( (uLW[0]+uLW[N])*(J[0]-J[N])-(uLW[N]+uLW[N-1])*(J[N]-J[N-1]) );



if (n==4500)  for(i=0;i<=N;i++){ file4<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.

if (n==5000)  for(i=0;i<=N;i++){ file5<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.

if (n==6000) for(i=0;i<=N;i++){ file6<<(-0.5*L+i*h)<<"  "<<uLWnew[i]<<endl; } //storing in the data file.



for(i=0;i<=N;i++){ u[i]=unew[i]; uL[i]=uLnew[i];  uLW[i]=uLWnew[i];  }


}


cout<<"end." <<endl; getch();

}